/**
 * @author Burcak Otlu Saritas
 *
 * 
 */

package calculations;

import java.io.BufferedReader;
import java.io.BufferedWriter;
import java.io.FileNotFoundException;
import java.io.FileReader;
import java.io.FileWriter;
import java.io.IOException;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import java.util.Map;
import java.util.Set;

import org.apache.commons.math3.distribution.BinomialDistribution;
import org.apache.commons.math3.util.ArithmeticUtils;
import org.apache.commons.math3.util.FastMath;

import annotate.encode.Commons;

public class CalculateBinomialDistributions {

	
	
	public void fillHashMapwithOccurences(Map<String,Integer> hashMap, String inputFileName){
		String strLine;
		FileReader fileReader = null;
		BufferedReader bufferedReader = null;
		
		int indexofFirstTab = 0;
		String name;
		Integer occurrence;
		
		try {
			fileReader = new FileReader(inputFileName);			
			bufferedReader = new BufferedReader(fileReader);
			
			while((strLine = bufferedReader.readLine())!=null) {
				
				indexofFirstTab = strLine.indexOf('\t');
				
				name = strLine.substring(0, indexofFirstTab);
				occurrence = Integer.valueOf(strLine.substring(indexofFirstTab+1));
				
				if (!(hashMap.containsKey(name))){
//					hashMap.put(strLine.toUpperCase(Locale.ENGLISH), Commons.LONG_ZERO);
					hashMap.put(name, occurrence);
				}
				
				strLine = null;
			}
			
			
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		} catch (IOException e) {
			e.printStackTrace();
		}
		
		try {
			bufferedReader.close();
			fileReader.close();
		} catch (IOException e) {
			e.printStackTrace();
		}			
	}
	
	public Long calculateWholeGenomeSize(List<Integer> chromosomeSizes, Long genomeSize){
		for(int i = 0; i<chromosomeSizes.size(); i++){
			genomeSize = genomeSize + chromosomeSizes.get(i);
		}
		return genomeSize;
	}
	
	public void initializeChromosomeSizes(List<Integer> chromosomeSizes){
		for(int i =0; i<Commons.NUMBER_OF_CHROMOSOMES_HG19; i++){
			chromosomeSizes.add(i, 0);
		}
	}
	
	public void getHg19ChromosomeSizes(List<Integer> chromosomeSizes, String inputFileName){
		String strLine = null;
		String chromName = null;
		int indexofFirstTab = 0;
		int chromSize = 0;
		
		FileReader fileReader;
		BufferedReader bufferedReader;
		
		try {
			fileReader = new FileReader(inputFileName);
			bufferedReader = new BufferedReader(fileReader);
			
			while((strLine = bufferedReader.readLine())!= null){
				indexofFirstTab = strLine.indexOf('\t');
				
				if (indexofFirstTab >= 0){
					chromName = strLine.substring(0, indexofFirstTab);
					chromSize = Integer.valueOf(strLine.substring(indexofFirstTab+1));
				}
				
				if (Commons.CHROMOSOME1.equals(chromName)){
					chromosomeSizes.set(0, chromSize);
				}else if (Commons.CHROMOSOME2.equals(chromName)){
					chromosomeSizes.set(1, chromSize);
				}else if (Commons.CHROMOSOME3.equals(chromName)){
					chromosomeSizes.set(2, chromSize);
				}else if (Commons.CHROMOSOME4.equals(chromName)){
					chromosomeSizes.set(3, chromSize);
				}else if (Commons.CHROMOSOME5.equals(chromName)){
					chromosomeSizes.set(4, chromSize);
				}else if (Commons.CHROMOSOME6.equals(chromName)){
					chromosomeSizes.set(5, chromSize);
				}else if (Commons.CHROMOSOME7.equals(chromName)){
					chromosomeSizes.set(6, chromSize);
				}else if (Commons.CHROMOSOME8.equals(chromName)){
					chromosomeSizes.set(7, chromSize);
				}else if (Commons.CHROMOSOME9.equals(chromName)){
					chromosomeSizes.set(8, chromSize);
				}else if (Commons.CHROMOSOME10.equals(chromName)){
					chromosomeSizes.set(9, chromSize);
				}else if (Commons.CHROMOSOME11.equals(chromName)){
					chromosomeSizes.set(10, chromSize);
				}else if (Commons.CHROMOSOME12.equals(chromName)){
					chromosomeSizes.set(11, chromSize);
				}else if (Commons.CHROMOSOME13.equals(chromName)){
					chromosomeSizes.set(12, chromSize);
				}else if (Commons.CHROMOSOME14.equals(chromName)){
					chromosomeSizes.set(13, chromSize);
				}else if (Commons.CHROMOSOME15.equals(chromName)){
					chromosomeSizes.set(14, chromSize);
				}else if (Commons.CHROMOSOME16.equals(chromName)){
					chromosomeSizes.set(15, chromSize);
				}else if (Commons.CHROMOSOME17.equals(chromName)){
					chromosomeSizes.set(16, chromSize);
				}else if (Commons.CHROMOSOME18.equals(chromName)){
					chromosomeSizes.set(17, chromSize);
				}else if (Commons.CHROMOSOME19.equals(chromName)){
					chromosomeSizes.set(18, chromSize);
				}else if (Commons.CHROMOSOME20.equals(chromName)){
					chromosomeSizes.set(19, chromSize);
				}else if (Commons.CHROMOSOME21.equals(chromName)){
					chromosomeSizes.set(20, chromSize);
				}else if (Commons.CHROMOSOME22.equals(chromName)){
					chromosomeSizes.set(21, chromSize);
				}else if (Commons.CHROMOSOMEX.equals(chromName)){
					chromosomeSizes.set(22, chromSize);
				}else if (Commons.CHROMOSOMEY.equals(chromName)){
					chromosomeSizes.set(23, chromSize);
				}				
			} // End of While
			
			bufferedReader.close();
			
		} catch (FileNotFoundException e) {
			e.printStackTrace();
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
	
	
	
	
	public void fillFunctionalElementList(Map<String,Integer> funcElementOccurrencesSearchInputHashMap,Map<String,Integer> funcElementOccurrencesWholeGenomeHashMap,Integer searchInputDataSize,long genomeSize, List<FunctionalElement> functionalElementList){
		
		Set<Map.Entry<String,Integer>> setSearchInput = funcElementOccurrencesSearchInputHashMap.entrySet();
		Iterator<Map.Entry<String, Integer>> itr = setSearchInput.iterator();
		double probability;
		
		while(itr.hasNext()){
			
			Map.Entry<String, Integer> entry = (Map.Entry<String,Integer>) itr.next();
			String functionalElementName = entry.getKey();
			Integer numberofOccurrencesinSearchInput = entry.getValue();
			
			//Get the corresponding numberofOccurrences of this functional Element from funcElementOccurrencesWholeGenomeHashMap
			
			int numberofOccurrencesinWholeGenome = funcElementOccurrencesWholeGenomeHashMap.get(functionalElementName);
			
			probability = (numberofOccurrencesinWholeGenome*1.0)/genomeSize;
				
			FunctionalElement functionalElement = new FunctionalElement(numberofOccurrencesinSearchInput,numberofOccurrencesinWholeGenome,probability,genomeSize,searchInputDataSize,functionalElementName);
			
			functionalElementList.add(functionalElement);
		}
		
	}
	
	
	
	public double logCombination(int n, int k){	
 		return  ArithmeticUtils.factorialLog(n)-ArithmeticUtils.factorialLog(k) -ArithmeticUtils.factorialLog(n-k);
	
	}
	

	public double func (int k, int n, double p){
		
		return logCombination(n,k) + (k *FastMath.log(p)) + ((n-k) *FastMath.log(1-p));
	}
	
	
	public double computePvalue(FunctionalElement element){
		double probability = element.getProbability();
		int  numberofOccurrencesinSearchInputData = element.getNumberofOccurrencesinSearchInputData();
		int searchInputSize = element.getSearchInputSize();
		double pValue = 0.0;
		
//		if (element.getName().equals("H3K27AC")){
//			System.out.println("here");
//		}
		
		for(int k = numberofOccurrencesinSearchInputData ; k<= searchInputSize; k++){
			pValue = pValue +  FastMath.pow(FastMath.E, func(k,searchInputSize,probability));
		}
				
		return pValue;
		
	}
	
	public double calculatePvalueUsingBinomialDistribution(FunctionalElement element){
		double probability = element.getProbability();
		int  numberofOccurrencesinSearchInputData = element.getNumberofOccurrencesinSearchInputData();
		int searchInputSize = element.getSearchInputSize();
		double pValue = 0.0;

		 BinomialDistribution binomialDistribution = new BinomialDistribution(searchInputSize, probability);
		 pValue = 1- binomialDistribution.cumulativeProbability(numberofOccurrencesinSearchInputData-1);
					
		return pValue;
	}
	
	public void fillPvalues(List<FunctionalElement> functionalElementList){
		for (int i = 0; i< functionalElementList.size(); i++){
			FunctionalElement element = functionalElementList.get(i);
		
			double pValue = computePvalue(element);
			element.setPValue(pValue);
			if (pValue<0.00001){
//				System.out.println(element.getName() + "\tp-value:\t" + pValue);
			}
		}
		
	}
	
	
	public void calculatePValuesUsingBinomialDistribution(List<FunctionalElement> functionalElementList){
		for (int i = 0; i< functionalElementList.size(); i++){
			FunctionalElement element = functionalElementList.get(i);
		
			double pValue = calculatePvalueUsingBinomialDistribution(element);
			element.setPValue(pValue);
		}
		
	}
	
	public void writePValues(List<FunctionalElement> functionalElementList, String outputFileName){
		FileWriter fileWriter = null;
		BufferedWriter bufferedWriter = null;
		
		try {
			fileWriter = new FileWriter(outputFileName);
			bufferedWriter = new BufferedWriter(fileWriter);
			
			for(int i = 0; i< functionalElementList.size(); i++){
				bufferedWriter.write(functionalElementList.get(i).getName() + "\t" + functionalElementList.get(i).getPValue()+ "\n");
				
			}
			
			bufferedWriter.close();
			
		} catch (IOException e) {
			e.printStackTrace();
		}
		
		
	}
	
	public void writeAll(List<FunctionalElement> dnaseCellLineNameList, String outputFileName){
		FileWriter fileWriter = null;
		BufferedWriter bufferedWriter = null;
		
		try {
			fileWriter = new FileWriter(outputFileName);
			bufferedWriter = new BufferedWriter(fileWriter);
			FunctionalElement element;

			//Header Row
			bufferedWriter.write("Name" + "\t"   + "OccurrencesInSearchInputData" + "\t" + "SearchInputDataSize" + "\t" + "OccurencesInWholeGenome" + "\t"+ "WholeGenomeSize" + "\t" + "PValue" + "\n");				

			for(int i = 0; i< dnaseCellLineNameList.size(); i++){
				element = dnaseCellLineNameList.get(i);
				bufferedWriter.write(element.getName() + "\t"   + element.getNumberofOccurrencesinSearchInputData() + "\t" +element.getSearchInputSize() + "\t" +element.getNumberofOccurrencesinWholeGenome() + "\t"+element.getGenomeSize() + "\t" + element.getPValue() + "\n");				
			}
			
			bufferedWriter.close();
			
		} catch (IOException e) {
			e.printStackTrace();
		}
	}
	
	
	
	public static void main(String[] args) {
		
		CalculateBinomialDistributions calculate = new CalculateBinomialDistributions();
		
		List<Integer> chromosomeSizes = new ArrayList<Integer>(24);
		Long genomeSize = Commons.LONG_ZERO;
		Integer searchInputDataSize = new Integer(2031);
		
		Map<String,Integer> dnaseCellLineNumberofNonoverlappingOccurrencesSearchInputHashMap = new HashMap<String,Integer>();
		Map<String,Integer> dnaseCellLineNumberofNonoverlappingOccurrencesWholeGenomeHashMap = new HashMap<String,Integer>();
		List<FunctionalElement> dnaseCellLineList = new ArrayList<FunctionalElement>();
		
		Map<String,Integer> tfbsNumberofNonoverlappingOccurrencesSearchInputHashMap = new HashMap<String,Integer>();
		Map<String,Integer> tfbsNumberofNonoverlappingOccurrencesWholeGenomeHashMap = new HashMap<String,Integer>();
		List<FunctionalElement> tfbsList = new ArrayList<FunctionalElement>();
		
		Map<String,Integer> histoneNumberofNonoverlappingOccurrencesSearchInputHashMap = new HashMap<String,Integer>();
		Map<String,Integer> histoneNumberofNonoverlappingOccurrencesWholeGenomeHashMap = new HashMap<String,Integer>();
		List<FunctionalElement> histoneList = new ArrayList<FunctionalElement>();
		
		//Calculate whole genome size using hg19_chromosome_sizes input file
		calculate.initializeChromosomeSizes(chromosomeSizes);
		calculate.getHg19ChromosomeSizes(chromosomeSizes,Commons.HG19_CHROMOSOME_SIZES_INPUT_FILE);
		genomeSize = calculate.calculateWholeGenomeSize(chromosomeSizes, genomeSize);
		
		
		//Dnase		
		calculate.fillHashMapwithOccurences(dnaseCellLineNumberofNonoverlappingOccurrencesSearchInputHashMap, Commons.NUMBER_OF_OCCURRENCES_OF_EACH_DNASE_CELL_LINE_SEARCH_INPUT_FILE);
		calculate.fillHashMapwithOccurences(dnaseCellLineNumberofNonoverlappingOccurrencesWholeGenomeHashMap, Commons.NUMBER_OF_OCCURRENCES_OF_EACH_DNASE_CELL_LINE_WHOLE_GENOME_OUTPUT_FILE);
		
		calculate.fillFunctionalElementList(dnaseCellLineNumberofNonoverlappingOccurrencesSearchInputHashMap,dnaseCellLineNumberofNonoverlappingOccurrencesWholeGenomeHashMap,searchInputDataSize,genomeSize,dnaseCellLineList);
		calculate.fillPvalues(dnaseCellLineList);
//		calculate.calculatePValuesUsingBinomialDistribution(dnaseCellLineList);
		calculate.writePValues(dnaseCellLineList, Commons.DNASE_CELLLINE_NAMES_P_VALUES);
		calculate.writeAll(dnaseCellLineList, Commons.DNASE_CELLLINE_NAMES_ALL_VALUES);
		
		//Tfbs
		calculate.fillHashMapwithOccurences(tfbsNumberofNonoverlappingOccurrencesSearchInputHashMap, Commons.NUMBER_OF_OCCURRENCES_OF_EACH_TFBS_SEARCH_INPUT_FILE);
		calculate.fillHashMapwithOccurences(tfbsNumberofNonoverlappingOccurrencesWholeGenomeHashMap, Commons.NUMBER_OF_OCCURRENCES_OF_EACH_TFBS_WHOLE_GENOME_OUTPUT_FILE);
		
		calculate.fillFunctionalElementList(tfbsNumberofNonoverlappingOccurrencesSearchInputHashMap,tfbsNumberofNonoverlappingOccurrencesWholeGenomeHashMap,searchInputDataSize,genomeSize,tfbsList);
		calculate.fillPvalues(tfbsList);
//		calculate.calculatePValuesUsingBinomialDistribution(tfbsList);

		calculate.writePValues(tfbsList, Commons.TFBS_P_VALUES);
		calculate.writeAll(tfbsList, Commons.TFBS_ALL_VALUES);
		
		//histone
		calculate.fillHashMapwithOccurences(histoneNumberofNonoverlappingOccurrencesSearchInputHashMap, Commons.NUMBER_OF_OCCURRENCES_OF_EACH_HISTONE_SEARCH_INPUT_FILE);
		calculate.fillHashMapwithOccurences(histoneNumberofNonoverlappingOccurrencesWholeGenomeHashMap, Commons.NUMBER_OF_OCCURRENCES_OF_EACH_HISTONE_WHOLE_GENOME_OUTPUT_FILE);
		
		calculate.fillFunctionalElementList(histoneNumberofNonoverlappingOccurrencesSearchInputHashMap,histoneNumberofNonoverlappingOccurrencesWholeGenomeHashMap,searchInputDataSize,genomeSize,histoneList);
		calculate.fillPvalues(histoneList);
//		calculate.calculatePValuesUsingBinomialDistribution(histoneList);

		calculate.writePValues(histoneList, Commons.HISTONE_P_VALUES);
		calculate.writeAll(histoneList, Commons.HISTONE_ALL_VALUES);
		
			
	}

}
